\subsection{重要抽样法}

\(I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) 中积分区域 \(C\) 可能是任意形状的,也可能无界, \(h\left( \mathbf{x}\right)\) 在 \(C\) 内各处的取值大小差异可能很大,使得直接用平均值法估计 \(I\) 时,很多样本点处于 \(\left| {h\left( \mathbf{x}\right) }\right|\) 接近于零的地方, 造成浪费,另外使得 \({\widehat{I}}_{2}\) 的渐近方差 (见(3.13)) 中的 \(\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right)\) 很大 \(\left( {\mathbf{U} \sim  \mathrm{U}\left( C\right) }\right)\) 。

{\color{red}\textbf{[类比]:} 重要抽样就像``有的放矢''——如果被积函数只在山峰处有值、平原处为零,均匀撒点就是浪费弹药。聪明的做法是在山峰密集投点、平原稀疏投点,然后用权重补偿这种不均匀性。}

为此,考虑用非均匀抽样: \(\left| {h\left( \mathbf{x}\right) }\right|\) 大的地方密集投点, \(\left| {h\left( \mathbf{x}\right) }\right|\) 小的地方稀疏投点。这样可以有效利用样本。设 \(g\left( \mathbf{x}\right) ,\mathbf{x} \in  C\) 是一个密度,其形状与 \(\left| {h\left( \mathbf{x}\right) }\right|\) 相近,且当 \(g\left( \mathbf{x}\right)  = 0\) 时 \(h\left( \mathbf{x}\right)  = 0\) ,当 \(\parallel \mathbf{x}\parallel  \rightarrow  \infty\) 时 \(h\left( \mathbf{x}\right)  = o\left( {g\left( \mathbf{x}\right) }\right)\) 。称 \(g\left( \mathbf{x}\right)\) 为试投密度或重要抽样密度。

设 \({\mathbf{X}}_{i}\) iid \(\sim  g\left( \mathbf{x}\right) ,i = 1,2,\ldots ,N\) 。令

\[
{\eta }_{i} = \frac{h\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) },i = 1,2,\ldots ,N
\]

则

\[
E{\eta }_{1} = {\int }_{C}\frac{h\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) }g\left( \mathbf{x}\right) d\mathbf{x} = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x} = I
\]

{\color{red}\textbf{[数学洞察]:} 重要抽样的核心公式$E[h(X)/g(X)] = \int h$看似魔术,实则是测度变换的直接应用——从$g$的测度下计算$h/g$的期望,等价于从Lebesgue测度下计算$h$的积分。}

因此可以用 \(\left\{  {{\eta }_{i},i = 1,2,\ldots ,N}\right\}\) 的样本平均值来估计 \(I\) ,即

\[
{\widehat{I}}_{3} = \bar{\eta } = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{h\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }. \tag{3.14}
\]

\({\widehat{I}}_{3}\) 是 \(I\) 的无偏估计和强相合估计。 \({\widehat{I}}_{3}\) 的渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \operatorname{Var}\left( \frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }\right) \frac{1}{N} = O\left( \frac{1}{N}\right) , \tag{3.15}
\]

{\color{red}\textbf{[方差推导]:} 这个等式基于蒙特卡洛估计的方差性质。设$\eta_i = h(X_i)/g(X_i)$为独立同分布的随机变量,则样本均值$\widehat{I}_3 = \frac{1}{N}\sum_{i=1}^N \eta_i$的方差为:
\[
\operatorname{Var}(\widehat{I}_3) = \operatorname{Var}\left(\frac{1}{N}\sum_{i=1}^N \eta_i\right) = \frac{1}{N^2} \cdot N \cdot \operatorname{Var}(\eta_i) = \frac{1}{N}\operatorname{Var}\left(\frac{h(\mathbf{X})}{g(\mathbf{X})}\right).
\]
由于$\operatorname{Var}(h(\mathbf{X})/g(\mathbf{X}))$是不依赖于$N$的常数(假设有限),因此$\operatorname{Var}(\widehat{I}_3) = O(1/N)$,表示方差以$1/N$的速率递减。这是蒙特卡洛方法的标准收敛速率。}

当 \(g\left( \mathbf{x}\right)\) 和 \(\left| {h\left( \mathbf{x}\right) }\right|\) 形状很接近时 \(\frac{\left| h\left( \mathbf{x}\right) \right| }{g\left( \mathbf{x}\right) }\) 近似为常数,方差 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 很小,这时 \({\widehat{I}}_{3}\) 的随机误差可以很小。

{\color{red}\textbf{[理想情况]:} 最优试投密度是$g(\mathbf{x}) \propto |h(\mathbf{x})|$,此时$h/g$确实为常数,方差为零!当然这需要知道$\int |h|$,相当于已经解决了问题。实际中我们只能用近似的$g$,但仍能大幅降低方差。}

用(3.14)估计 \(I\) 与 \(§{2.2.4}\) 的舍选法 II 有类似的想法,这种方法叫做重要抽样法 (importance sampling), 是随机模拟的重要方法。

为什么当 \(g\left( \mathbf{x}\right)\) 和 \(\left| {h\left( \mathbf{x}\right) }\right|\) 形状很接近时 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 很小？事实上,

{\color{red}\textbf{[证明目标]:} 我们要严格证明:方差$\operatorname{Var}(\widehat{I}_3)$有下界,且当$g(\mathbf{x}) \propto |h(\mathbf{x})|$时达到最小值。证明策略是利用Jensen不等式建立方差的下界。}

\[
\operatorname{Var}\left( \frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }\right)  = E\left( \frac{{h}^{2}\left( \mathbf{X}\right) }{{g}^{2}\left( \mathbf{X}\right) }\right)  - \left( {E\frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }}\right)  = E\left( \frac{{h}^{2}\left( \mathbf{X}\right) }{{g}^{2}\left( \mathbf{X}\right) }\right)  - {\left( {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\right) }^{2}
\]

{\color{red}\textbf{[为什么用Jensen不等式]:} 注意到$E(h^2/g^2)$很难直接处理,但我们可以利用$h^2 = (|h|)^2$,然后对凸函数$\phi(t) = t^2$应用Jensen不等式。Jensen不等式说:若$\phi$是凸函数,则$E[\phi(Y)] \geq \phi(E[Y])$。这里取$Y = |h(\mathbf{X})|/g(\mathbf{X})$,则$\phi(Y) = Y^2 = h^2/g^2$,于是$E(h^2/g^2) \geq (E|h|/g)^2$。}

{\color{red}\textbf{[数学推导第一步]:} 应用Jensen不等式:由于$t \mapsto t^2$是凸函数,对随机变量$|h(\mathbf{X})|/g(\mathbf{X})$有
\[E\left(\frac{h^2(\mathbf{X})}{g^2(\mathbf{X})}\right) = E\left[\left(\frac{|h(\mathbf{X})|}{g(\mathbf{X})}\right)^2\right] \geq \left[E\left(\frac{|h(\mathbf{X})|}{g(\mathbf{X})}\right)\right]^2.\]
因此方差$\geq (E|h|/g)^2 - I^2$。}

\[
\geq  {\left( E\frac{\left| h\left( \mathbf{X}\right) \right| }{g\left( \mathbf{X}\right) }\right) }^{2} - {I}^{2}\;\text{ (由 Jensen 不等式) } \tag{3.16}
\]

{\color{red}\textbf{[数学推导第二步]:} 计算$E(|h(\mathbf{X})|/g(\mathbf{X}))$:注意$\mathbf{X} \sim g(\mathbf{x})$,所以
\[E\left(\frac{|h(\mathbf{X})|}{g(\mathbf{X})}\right) = \int_C \frac{|h(\mathbf{x})|}{g(\mathbf{x})} \cdot g(\mathbf{x}) d\mathbf{x} = \int_C |h(\mathbf{x})| d\mathbf{x}.\]
这就是被积函数绝对值的积分!}

\[
= {\left( {\int }_{C}\left| h\left( \mathbf{x}\right) \right| d\mathbf{x}\right) }^{2} - {I}^{2}, \tag{3.17}
\]

{\color{red}\textbf{[最优条件分析]:} Jensen不等式取等号的条件是随机变量$|h(\mathbf{X})|/g(\mathbf{X})$为常数(几乎处处)。设这个常数为$k$,即$|h(\mathbf{x})| = k \cdot g(\mathbf{x})$对所有$\mathbf{x} \in C$成立。对两边积分:$\int_C |h(\mathbf{x})| d\mathbf{x} = k \int_C g(\mathbf{x}) d\mathbf{x} = k$。因此$k = \int_C |h(\mathbf{x})| d\mathbf{x}$,所以最优试投密度为$g(\mathbf{x}) = |h(\mathbf{x})|/k = |h(\mathbf{x})|/\int_C |h(\mathbf{t})| dt$。}

{\color{red}\textbf{[直觉理解]:} 为什么$g \propto |h|$最优?因为此时每个样本点$\mathbf{X}_i$都``等价重要'':虽然$h(\mathbf{X}_i)/g(\mathbf{X}_i)$的绝对值是常数,但符号可能不同(如果$h$有正有负)。方差来自$h$的符号变化,而不是幅度变化。如果$h \geq 0$,则$h/g$完全为常数,方差为零!}

当且仅当 \(\frac{\left| h\left( \mathbf{X}\right) \right| }{g\left( \mathbf{X}\right) }\) 为常数时不等式(3.16)中的等号成立,即 \(g\left( \mathbf{x}\right)  = \frac{1}{{\int }_{C}\left| {h\left( \mathbf{t}\right) }\right| {dt}}\left| {h\left( \mathbf{x}\right) }\right| ,\mathbf{x} \in  C\) 时 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 达到最小。

上述求积分的问题也可以表述为求随机变量函数的期望问题。设 \(\mathbf{Y} \sim  f\left( \mathbf{y}\right)\) ,要求 \(\mathbf{Y}\) 的函数 \(h\left( \mathbf{Y}\right)\) 的期望值 \({Eh}\left( \mathbf{Y}\right)  = \int h\left( \mathbf{y}\right) f\left( \mathbf{y}\right) d\mathbf{y}\) ,可以抽取 \(\mathbf{Y}\) 的随机数 \({\mathbf{Y}}_{i},i = 1,2,\ldots ,N\) ,然后用平均值法 \(\frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {\mathbf{Y}}_{i}\right)\) 估计 \({Eh}\left( \mathbf{Y}\right)\) 。但是,如果很难生成 \(\mathbf{Y}\) 的随机数,或者 \(\mathbf{Y}\) 的分布集中于 \(h\left( \mathbf{y}\right)\) 接近于零的位置以至于积分效率很低,可以找一个试投密度 \(g\left( \mathbf{x}\right)\) ,从 \(g\left( \mathbf{x}\right)\) 产生随机数 \({\mathbf{X}}_{i},i = 1,2,\ldots ,N\) ,用如下重要抽样法

\[
{\widehat{I}}_{3.1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {\mathbf{X}}_{i}\right) \frac{f\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) } \tag{3.18}
\]

估计 \({Eh}\left( \mathbf{Y}\right)\) ,易见 \(E{\widehat{I}}_{3.1} = {Eh}\left( \mathbf{Y}\right)\) 。称 \({W}_{i} = \frac{f\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) 为重要性权重, \({Eh}\left( \mathbf{Y}\right)\) 可以用重要性权重估计为

\[
{\widehat{I}}_{3.1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) . \tag{3.19}
\]

选取试投密度 \(g\left( \mathbf{x}\right)\) 时,要求当 \(h\left( \mathbf{x}\right) f\left( \mathbf{x}\right)  \neq  0\) 时 \(g\left( \mathbf{x}\right)  \neq  0\) ,当 \(\mathbf{x}\) 趋于无穷时 \(h\left( \mathbf{x}\right) f\left( \mathbf{x}\right)  =\)  \(o\left( {g\left( \mathbf{x}\right) }\right)\) ( \(g\left( \mathbf{x}\right)\) 相对厚尾),一般还要求 \(g\left( \mathbf{x}\right)\) 的形状与 \(\left| {h\left( \mathbf{x}\right) }\right| f\left( \mathbf{x}\right)\) 形状接近。如果 \(g\left( \mathbf{x}\right)\) 的相对厚尾性难以确定, 可以使用如下保险的试投密度

\[
\widetilde{g}\left( \mathbf{x}\right)  = {\rho g}\left( \mathbf{x}\right)  + \left( {1 - \rho }\right) r\left( \mathbf{x}\right) , \tag{3.20}
\]

其中 \(\rho\) 接近于 1, \(r\left( \mathbf{x}\right)\) 是柯西分布或 Pareto 分布这样的重尾分布。要生成 \(N\) 个 \(\widetilde{g}\left( \mathbf{x}\right)\) 的随机数,只要生成 \({N\rho }\) 个 \(g\left( \mathbf{x}\right)\) 的随机数和 \(N\left( {1 - \rho }\right)\) 个 \(r\left( \mathbf{x}\right)\) 的随机数。

选取不适当的试投密度会把绝大多数样本点投到了对计算积分不重要的位置, 使得样本点中只有极少数点是真正有作用的。在多维问题中合适的试投密度尤其难找, 经常需要反复试验。

有时 \(\mathbf{Y} \sim  f\left( \cdot \right)\) 不仅很难直接抽样,而且 \(f\left( \cdot \right)\) 本身未知,只能确定到差一个常数倍的 \(\widetilde{f}\left( \mathbf{x}\right)  = {cf}\left( \mathbf{x}\right)\) ,常数 \(c\) 未知,为了求常数 \(c\) 需要计算 \(c = \int \widetilde{f}\left( \mathbf{y}\right) d\mathbf{y}\) ,计算 \(c\) 一般很困难。这时,定义重要性权重为 \({W}_{i} = \frac{\widetilde{f}\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) ,公式(3.19)可以改成

\[
{\widehat{I}}_{4} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}} \tag{3.21}
\]

这称为标准化重要抽样法。对(3.21)的分子和分母都除以 \({cN}\) 后分子 a.s. 收敛到 \({Eh}\left( \mathbf{Y}\right)\) ,分母 a.s. 收敛到 1,所以 (3.21) 是 \(\operatorname{Eh}\left( \mathbf{Y}\right)\) 的强相合估计,但不是无偏的。标准化重要抽样估计往往比无偏估计 \({\widehat{I}}_{3.1}\) 有更小的均方误差。关于标准化重要抽样法的渐近方差的讨论参见 \(\operatorname{Liu}{\left( {2001}\right) }^{\left\lbrack  {28}\right\rbrack  }§{2.5.3}\) 。

如果需要对多个不同的函数 \(h\left( \cdot \right)\) 计算 \({Eh}\left( \mathbf{Y}\right)\) ,则选取试抽样密度 \(g\left( \mathbf{x}\right)\) 时应使得 \(g\left( \mathbf{x}\right)\) 尽可能与 \(\mathbf{Y}\) 的密度 \(f\left( \mathbf{x}\right)\) 形状接近,这样权重 \({W}_{i} = \frac{\widehat{f}\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) 的分布不过于偏斜,不至于出现绝大部分权重集中于少数样本点的情形。抽样值 \({\mathbf{X}}_{i}\) 与权重 \({W}_{i}\) 一起可以看作是分布 \(f\left( \cdot \right)\) 的某种抽样。

{\color{red}\textbf{[背景说明]:} 前面的重要抽样中,样本点$\mathbf{X}_i$来自试投密度$g$,通过权重$W_i = f(\mathbf{X}_i)/g(\mathbf{X}_i)$来``纠正''分布。但实际中我们可能用更复杂的方法产生$({\mathbf{X}_i, W_i})$对(如MCMC、粒子滤波)。``适当加权抽样''定义了什么样的加权样本是``有效的''。}

定义 (适当加权抽样) 随机变量序列 \(\left\{  {\left( {{\mathbf{X}}_{i},{W}_{i}}\right) ,i = 1,2,\ldots ,N}\right\}\) 称为关于密度 \(f\left( \cdot \right)\) 的适当加权抽样,如果对于任何平方可积函数 \(h\left( \cdot \right)\) 都有

\[
E\left\lbrack  {h\left( {X}_{i}\right) {W}_{i}}\right\rbrack   = c{E}_{f}\left\lbrack  {h\left( X\right) }\right\rbrack   = c\int h\left( \mathbf{x}\right) f\left( \mathbf{x}\right) d\mathbf{x},i = 1,2,\ldots ,N,
\]

其中 \(c\) 是归一化常数。

{\color{red}\textbf{[定义的直觉理解]:} ``适当加权''的核心思想:虽然样本点$\mathbf{X}_i$不是直接从目标密度$f$抽取的,但通过权重$W_i$加权后,期望$E[h(\mathbf{X}_i)W_i]$与目标期望$E_f[h(X)]$成正比(相差常数$c$)。这意味着用$\sum h(\mathbf{X}_i)W_i / \sum W_i$可以一致地估计$E_f[h(X)]$,对所有函数$h$都成立!}

{\color{red}\textbf{[为什么需要常数$c$]:} 常数$c$的存在是因为:1) 实际中$f$可能只知道到差一个常数($\widetilde{f} = cf$),2) 权重$W_i$的尺度可以任意缩放而不影响标准化估计$\sum h(\mathbf{X}_i)W_i / \sum W_i$。关键是比例关系,而非绝对值。}

设随机变量(X, W)联合密度为 \(g\left( {\mathbf{x},w}\right)\) ,则(X, W)的样本为密度 \(f\left( \cdot \right)\) 的适当加权抽样的充分必要条件是

\[
{E}_{g}\left( {W \mid  \mathbf{x}}\right)  = {E}_{g}\left( W\right) \frac{f\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) },\forall \mathbf{x},
\]

其中 \({E}_{g}\left( W\right)\) 是关于 \(W\) 的边缘密度的期望, \({E}_{g}\left( {W \mid  \mathbf{x}}\right)\) 是在(X, W)的联合密度下条件期望 \(E\left( {W \mid  \mathbf{X}}\right)\) 在 \(\mathbf{X} = \mathbf{x}\) 处的值。

{\color{red}\textbf{[充要条件的含义]:} 这个条件说:给定位置$\mathbf{x}$,权重$W$的条件期望必须正比于似然比$f(\mathbf{x})/g(\mathbf{x})$,比例系数是$E_g(W)$。直观地:在$\mathbf{x}$处,平均权重应该反映$f$相对于$g$的``超出程度''。}

{\color{red}\textbf{[为什么这个条件保证适当加权]:} 证明思路:计算$E[h(\mathbf{X})W]$。设$(\mathbf{X}, W)$的联合密度为$g(\mathbf{x}, w)$,则
\begin{align*}
E[h(\mathbf{X})W] &= \iint h(\mathbf{x}) w \cdot g(\mathbf{x}, w) d\mathbf{x} dw \\
&= \int h(\mathbf{x}) \left[\int w \cdot g(\mathbf{x}, w) dw\right] d\mathbf{x} \\
&= \int h(\mathbf{x}) E_g(W|\mathbf{x}) g(\mathbf{x}) d\mathbf{x} \quad \text{(这里$g(\mathbf{x})$是$\mathbf{X}$的边缘密度)}\\
&= \int h(\mathbf{x}) \cdot E_g(W) \frac{f(\mathbf{x})}{g(\mathbf{x})} \cdot g(\mathbf{x}) d\mathbf{x} \quad \text{(用充要条件)}\\
&= E_g(W) \int h(\mathbf{x}) f(\mathbf{x}) d\mathbf{x} = E_g(W) \cdot E_f[h(X)].
\end{align*}
取$c = E_g(W)$即得定义!反之,若对所有$h$都有$E[h(\mathbf{X})W] = c E_f[h(X)]$,反推可得充要条件。}

{\color{red}\textbf{[经典例子]:} 标准重要抽样:$\mathbf{X} \sim g(\mathbf{x})$,权重$W = f(\mathbf{X})/g(\mathbf{X})$是$\mathbf{X}$的函数(确定性)。此时$E_g(W|\mathbf{x}) = f(\mathbf{x})/g(\mathbf{x})$ (给定$\mathbf{x}$,$W$就确定了),而$E_g(W) = \int [f(\mathbf{x})/g(\mathbf{x})] g(\mathbf{x}) d\mathbf{x} = \int f(\mathbf{x}) d\mathbf{x} = 1$。所以$E_g(W|\mathbf{x}) = 1 \cdot f(\mathbf{x})/g(\mathbf{x}) = E_g(W) \cdot f(\mathbf{x})/g(\mathbf{x})$,满足充要条件!}

{\color{red}\textbf{[权重退化问题]:} 重要抽样的主要问题:当试投密度$g$与目标密度$f$差异较大时,权重$W_i = f(\mathbf{X}_i)/g(\mathbf{X}_i)$会高度不均匀。比如99\%的样本权重接近0,只有1\%的样本有显著权重。这导致``有效样本量''远小于名义样本量$N$,估计效率很低。这就是权重退化(weight degeneracy)现象。}

在重要抽样法和标准化重要抽样法的实际应用中, 好的试抽样分布很难获得, 所以权重 \(\left\{  {{W}_{i} = f\left( {\mathbf{X}}_{i}\right) /g\left( {\mathbf{X}}_{i}\right) }\right\}\) 经常会差别很大,使得抽样样本主要集中在少数几个权重最大的样本点上。为此, 可以舍弃权重太小的样本点, 重新抽样替换这样的样本点。这种方法称为舍选控制 (Rejection Control),描述如下。

{\color{red}\textbf{[舍选控制的基本思想]:} 核心策略:1) 对权重太小($W_i < c$)的样本点,以概率$W_i/c$接受,否则舍弃并重新抽取。这类似于舍选法(rejection sampling)!2) 对被接受的样本,调整其权重,使得$({\mathbf{X}_i, W_i^*})$仍是适当加权抽样。目标:减少权重的变异性,提高有效样本量。}

首先,需要选定权重的一个阈值 \(c\) ,然后对每个样本点 \({\mathbf{X}}_{i}\) ,若 \({W}_{i} \geq  c\) ,则接受 \({\mathbf{X}}_{i}\) ; 若 \({W}_{i} < c\) ,则以概率 \({W}_{i}/c\) 接受 \({\mathbf{X}}_{i}\) ,否则舍弃 \({\mathbf{X}}_{i}\) ,重新抽取。最后,所有样本点的权重调整为新的

{\color{red}\textbf{[算法流程详解]:} 对每个从$g$抽取的$\mathbf{X}_i$,计算$W_i = f(\mathbf{X}_i)/g(\mathbf{X}_i)$:
\begin{itemize}
\item 若$W_i \geq c$:必然接受。这些是``重要''样本,目标密度$f$在此处相对较大。
\item 若$W_i < c$:以概率$W_i/c < 1$接受。权重越小,接受概率越低。若不接受,舍弃此样本,重新从$g$抽取新样本,重复此流程。
\end{itemize}
这样,权重小的样本被``稀疏化'',权重大的样本被保留,从而减少权重的不平衡。}

\[
{W}_{i}^{ * } = {p}_{c}\frac{{W}_{i}}{\min \left\{  {1,{W}_{i}/c}\right\}  } = {p}_{c}\max \left\{  {{W}_{i},c}\right\}  , \tag{3.22}
\]

{\color{red}\textbf{[公式(3.22)解释]:} 调整后的权重$W_i^* = p_c \max\{W_i, c\}$的含义:
\begin{itemize}
\item 若$W_i \geq c$:$W_i^* = p_c W_i$(乘以归一化常数)。
\item 若$W_i < c$:$W_i^* = p_c c$。注意这个样本是以概率$W_i/c$被接受的,所以其权重需要``放大''为$c$(再乘$p_c$)来补偿接受概率的降低。
\end{itemize}
直觉:被接受的低权重样本代表了多个被舍弃的样本,所以权重需要上调。等价变换:$W_i / \min\{1, W_i/c\} = W_i / (W_i/c) = c$ (当$W_i < c$时)。}

{\color{red}\textbf{[为什么需要归一化常数$p_c$]:} 由于舍选过程改变了$\mathbf{X}_i$的分布(从$g$变为$g^*$),需要$p_c$来保证$({\mathbf{X}_i, W_i^*})$仍是适当加权抽样。$p_c$本质上是接受概率的期望,即从$g$抽取一个样本被接受的概率。}

其中

\[
{p}_{c} = \int \min \left\{  {1,\frac{f\left( \mathbf{x}\right) }{{cg}\left( \mathbf{x}\right) }}\right\}  g\left( \mathbf{x}\right) d\mathbf{x} = \int \min \left\{  {g\left( \mathbf{x}\right) ,\frac{f\left( \mathbf{x}\right) }{c}}\right\}  d\mathbf{x} \tag{3.23}
\]

{\color{red}\textbf{[公式(3.23)推导]:} $p_c$的计算:从$g$抽取$\mathbf{x}$,接受概率为
\[a(\mathbf{x}) = \begin{cases} 1, & \text{if } W(\mathbf{x}) = f(\mathbf{x})/g(\mathbf{x}) \geq c \\ W(\mathbf{x})/c = f(\mathbf{x})/(cg(\mathbf{x})), & \text{if } W(\mathbf{x}) < c \end{cases} = \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\}.\]
所以总接受概率$p_c = E_{g}[a(\mathbf{X})] = \int a(\mathbf{x}) g(\mathbf{x}) d\mathbf{x} = \int \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\} g(\mathbf{x}) d\mathbf{x} = \int \min\{g(\mathbf{x}), f(\mathbf{x})/c\} d\mathbf{x}$。}

是归一化常数,如果使用标准化重要抽样法, \({p}_{c}\) 可以省略,否则, \({p}_{c}\) 可以估计为

\[
{\widehat{p}}_{c} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\min \left\{  {1,\frac{{W}_{i}}{c}}\right\}  . \tag{3.24}
\]

这样得到的 \(\left\{  {\left( {{\mathbf{X}}_{i},{W}_{i}^{ * }}\right) ,i = 1,\ldots ,N}\right\}\) 是关于 \(f\left( \cdot \right)\) 适当加权的,被接受的 \({\mathbf{X}}_{i}\) 的分布密度为

\[
{g}^{ * }\left( \mathbf{x}\right)  = \frac{1}{{p}_{c}}\min \left\{  {g\left( \mathbf{x}\right) ,\frac{f\left( \mathbf{x}\right) }{c}}\right\}  . \tag{3.25}
\]

{\color{red}\textbf{[公式(3.25)推导]:} 接受密度$g^*(\mathbf{x})$的推导:从$g(\mathbf{x})$抽取$\mathbf{x}$,以概率$a(\mathbf{x}) = \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\}$接受。根据条件概率,被接受样本的密度为
\[g^*(\mathbf{x}) = \frac{g(\mathbf{x}) \cdot a(\mathbf{x})}{p_c} = \frac{g(\mathbf{x}) \cdot \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\}}{p_c} = \frac{\min\{g(\mathbf{x}), f(\mathbf{x})/c\}}{p_c}.\]
注意$g^*$自动满足$\int g^*(\mathbf{x}) d\mathbf{x} = 1$(因为$p_c$就是归一化常数)。}

{\color{red}\textbf{[验证适当加权性]:} 需要验证$({\mathbf{X}_i, W_i^*})$ (其中$\mathbf{X}_i \sim g^*$,$W_i^* = p_c \max\{W_i, c\}$)是适当加权抽样。计算:当$\mathbf{X} \sim g^*$时,
\begin{align*}
E[h(\mathbf{X})W^*] &= \int h(\mathbf{x}) \cdot p_c \max\{f(\mathbf{x})/g(\mathbf{x}), c\} \cdot g^*(\mathbf{x}) d\mathbf{x} \\
&= \int h(\mathbf{x}) \max\{f(\mathbf{x})/g(\mathbf{x}), c\} \cdot \min\{g(\mathbf{x}), f(\mathbf{x})/c\} d\mathbf{x}.
\end{align*}
分两种情况:1) 若$f(\mathbf{x})/g(\mathbf{x}) \geq c$,则$\max = f/g$,$\min = g$,乘积$= f$。2) 若$f(\mathbf{x})/g(\mathbf{x}) < c$,则$\max = c$,$\min = f/c$,乘积$= f$。因此$E[h(\mathbf{X})W^*] = \int h(\mathbf{x}) f(\mathbf{x}) d\mathbf{x} = E_f[h(X)]$,满足适当加权(常数$c=1$)!}

{\color{red}\textbf{[方法总结]:} 舍选控制通过``舍弃+权重调整''实现了两个目标:1) 减少权重的变异性(所有调整后权重$\geq p_c c$),2) 保持适当加权性质。代价是需要重新抽取被舍弃的样本,增加了计算量。阈值$c$的选择是关键:太小则舍弃太少,效果不明显;太大则舍弃太多,计算量大增。实践中常取权重的中位数或某个分位数。}

阈值 \(c\) 在实际中可以从权重 \(\left\{  {W}_{i}\right\}\) 选取,比如,取为 \(\left\{  {W}_{i}\right\}\) 的某个分位数。

例 3.2.3. 用 \(\mathrm{{MC}}\) 积分法计算 \(I = {\int }_{0}^{1}{e}^{x}{dx} = e - 1 \approx  {1.718}\) 。对被积函数 \(h\left( x\right)  = {e}^{x}\) 做泰勒展开得

\[
{e}^{x} = 1 + x + \frac{{x}^{2}}{2!} + \cdots
\]

取

\[
g\left( x\right)  = c\left( {1 + x}\right)  = \frac{2}{3}\left( {1 + x}\right)
\]

要产生 \(g\left( x\right)\) 的随机数可以用逆变换法,密度 \(g\left( x\right)\) 的分布函数 \(G\left( x\right)\) 的反函数为

\[
{G}^{-1}\left( y\right)  = \sqrt{1 + {3y}} - 1,0 < y < 1
\]

因此,取 \({U}_{i}\) iid \(\mathrm{U}\left( {0,1}\right)\) ,令 \({X}_{i} = \sqrt{1 + 3{U}_{i}} - 1,i = 1,2,\ldots ,N\) ,则重要抽样法的积分公式为

\[
{\widehat{I}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{{e}^{{X}_{i}}}{\frac{2}{3}\left( {1 + {X}_{i}}\right) }
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \frac{1}{N}\left( {\frac{3}{2}{\int }_{0}^{1}\frac{{e}^{2x}}{1 + x}{dx} - {I}^{2}}\right)  \approx  {0.02691}/N.
\]

如果用平均值法, 估计公式为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{e}^{{U}_{i}}
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{1}{N}\left( {{\int }_{0}^{1}{e}^{2x}{dx} - {I}^{2}}\right)  \approx  {0.2420}/N \approx  {9.0} \times  \operatorname{Var}\left( {\widehat{I}}_{3}\right)
\]

是重要抽样法方差的 9 倍。

如果用随机投点法, \(h\left( x\right)  = {e}^{x} \leq  e\left( {0 < x < 1}\right)\) ,取上界 \(M = e\) ,向 \(\left\lbrack  {0,1}\right\rbrack   \times  \left\lbrack  {0,M}\right\rbrack\) 随机投点,落到 \(f\left( x\right)\) 下方的概率为

\[
p = I/\left( {M\left( {b - a}\right) }\right)  = \left( {e - 1}\right) /e,
\]

设投 \(N\) 点落到 \(h\left( x\right)\) 下方的频率为 \(\widehat{p}\) ,用随机投点法估计 \(I\) 的公式为

\[
{\widehat{I}}_{1} = \widehat{p} \cdot  M\left( {b - a}\right)  = e\widehat{p},
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = {e}^{2}p\left( {1 - p}\right) /N = \left( {e - 1}\right) /N \approx  {1.718}/N \approx  {7.1} \times  \operatorname{Var}\left( {\widehat{I}}_{2}\right)  \approx  {64.8} \times  \operatorname{Var}\left( {\widehat{I}}_{3}\right)
\]

可见选择合适的抽样算法对减少计算量、提高精度是十分重要的。

例 3.2.4. 设二元函数 \(f\left( {x,y}\right)\) 定义如下

\[
f\left( {x,y}\right)  = \exp \left\{  {-{45}{\left( x + {0.4}\right) }^{2} - {60}{\left( y - {0.5}\right) }^{2}}\right\}   + {0.5}\exp \left\{  {-{90}{\left( x - {0.5}\right) }^{2} - {45}{\left( y + {0.1}\right) }^{4}}\right\}
\]

求如下二重定积分

\[
I = {\int }_{-1}^{1}{\int }_{-1}^{1}f\left( {x,y}\right) {dxdy}
\]

\(f\left( {x,y}\right)\) 有两个分别以(-0.4,0.5)和(0.5, - 0.1)为中心的峰,对积分有贡献的区域主要集中在 (-0.4,0.5)和(0.5, - 0.1)附近,在其他地方函数值很小,对积分贡献很小。

用平均值法 (3.12),取点数 \(N = {10000}\) 时 \({\widehat{I}}_{2}\) 的一个估计值为 \({\widehat{I}}_{2} = {0.132}\) ,从这一次模拟估计的 \({\widehat{I}}_{2}\) 的渐近方差为 \({1.86} \times  {10}^{-5}\) 。重复模拟 \(B = {100}\) 次,得到这 100 个 \({\widehat{I}}_{2}\) 的平均值为 0.1256,样本方差为 \({1.95} \times  {10}^{-5}\) 。

用重要抽样法。取试投密度为

\[
g\left( {x,y}\right)  \propto  \widetilde{g}\left( {x,y}\right)
\]

\[
= \exp \left\{  {-{45}{\left( x + {0.4}\right) }^{2} - {60}{\left( y - {0.5}\right) }^{2}}\right\}   + {0.5}\exp \left\{  {-{90}{\left( x - {0.5}\right) }^{2} - {10}{\left( y + {0.1}\right) }^{2}}\right\}  ,
\]

\[
- \infty  < x < \infty , - \infty  < y < \infty ,
\]

这样抽取到 \(\left\lbrack  {-1,1}\right\rbrack   \times  \left\lbrack  {-1,1}\right\rbrack\) 范围外的点对积分没有贡献,因为构成 \(g\left( {x,y}\right)\) 的两个密度都很集中,所以效率损失不大。需要求使得 \(\widetilde{g}\left( {x,y}\right)\) 化为密度的比例常数。记 \(\mathrm{N}\left( {\mu ,{\sigma }^{2}}\right)\) 的分布密度为 \(f\left( {x;\mu ,{\sigma }^{2}}\right)\) ,对 \(\widetilde{g}\left( {x,y}\right)\) 积分得

\[
{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\widetilde{g}\left( {x,y}\right)  = \sqrt{{2\pi }/{90}}{\int }_{-\infty }^{\infty }f\left( {x; - {0.4},{90}^{-1}}\right) {dx} \cdot  \sqrt{{2\pi }/{120}}{\int }_{-\infty }^{\infty }f\left( {y;{0.5},{120}^{-1}}\right) {dy}
\]

\[
+ {0.5}\sqrt{{2\pi }/{180}}{\int }_{-\infty }^{\infty }f\left( {x;{0.5},{180}^{-1}}\right) {dx} \cdot  \sqrt{{2\pi }/{20}}{\int }_{-\infty }^{\infty }f\left( {y; - {0.1},{20}^{-1}}\right) {dy}
\]

\[
= \sqrt{{2\pi }/{90}}\sqrt{{2\pi }/{120}} + {0.5}\sqrt{{2\pi }/{180}}\sqrt{{2\pi }/{20}}
\]

\[
\approx  {0.1128199}
\]

于是令

\[
g\left( {x,y}\right)  = \widetilde{g}\left( {x,y}\right) /{0.1128199}
\]

\[
= {0.5358984f}\left( {x; - {0.4},{90}^{-1}}\right) f\left( {y;{0.5},{120}^{-1}}\right)
\]

\[
+ {0.4641016f}\left( {x;{0.5},{180}^{-1}}\right) f\left( {y; - {0.1},{20}^{-1}}\right) \text{ , }
\]

\[
- \infty  < x < \infty , - \infty  < y < \infty ,
\]

用复合抽样法对 \(g\left( {x,y}\right)\) 抽样, \(N = {10000}\) 时一次模拟得到的 \({\widehat{I}}_{3} = {0.126}\) ,从这一次模拟估计的 \({\widehat{I}}_{3}\) 的渐近方差为 \({6.84} \times  {10}^{-8}\) ,重复模拟 \(B = {100}\) 次,则 100 个 \({\widehat{I}}_{3}\) 的平均值为 0.1258,样本方差为 \({5.37} \times  {10}^{-8}\) 。 \({\widehat{I}}_{2}\) 的样本方差是 \({\widehat{I}}_{3}\) 的样本方差的 363 倍,如果要达到相同的估计方差, 两种方法的样本量相差三百多倍。

例 3.2.5. 标准化的重要抽样法在贝叶斯统计推断中有重要作用。例如,设独立的观测样本 \({Y}_{j}\) 服从如下的贝塔一二项分布:

\[
f\left( {{y}_{j} \mid  K,\eta }\right)  = P\left( {{Y}_{j} = {y}_{j}}\right)  = \left( \begin{matrix} {n}_{j} \\  {y}_{j} \end{matrix}\right) \frac{B\left( {{K\eta } + {y}_{j},K\left( {1 - \eta }\right)  + {n}_{j} - {y}_{j}}\right) }{B\left( {{K\eta },K\left( {1 - \eta }\right) }\right) },{y}_{j} = 0,1,\ldots ,{n}_{j}, \tag{3.26}
\]

其中 \(B\left( {\cdot , \cdot  }\right)\) 是贝塔函数, \({n}_{j}\) 为已知的正整数, \(K > 0,0 < \eta  < 1\) 为未知参数。贝塔一二项分布用于描述比二项分布更为分散的随机变量分布。按照贝叶斯统计的做法,假设参数 \(\left( {K,\eta }\right)\) 也是随机变量,具有所谓的 “先验分布”,假设 \(\left( {K,\eta }\right)\) 有如下的 “无信息” 先验分布密度:

\[
\pi \left( {K,\eta }\right)  \propto  \frac{1}{{\left( 1 + K\right) }^{2}}\frac{1}{\eta \left( {1 - \eta }\right) }, \tag{3.27}
\]

则 \(\left( {K,\eta }\right)\) 有如下的 “后验密度”:

\[
\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right)  \propto  \pi \left( {K,\eta }\right) \mathop{\prod }\limits_{{j = 1}}^{n}f\left( {{y}_{j} \mid  K,\eta }\right)
\]

\[
\propto  \frac{1}{{\left( 1 + K\right) }^{2}}\frac{1}{\eta \left( {1 - \eta }\right) }\mathop{\prod }\limits_{{j = 1}}^{n}\frac{B\left( {{K\eta } + {y}_{j},K\left( {1 - \eta }\right)  + {n}_{j} - {y}_{j}}\right) }{B\left( {{K\eta },K\left( {1 - \eta }\right) }\right) }. \tag{3.28}
\]

设要求 \(E\left( {\log K \mid  \mathbf{Y}}\right)  = {\int }_{0}^{\infty }\log K\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right) {dK}\) 的值。

如果可以从后验密度 \(\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right)\) 直接抽样,可以用平均值法估计 \(E\left( {\log K \mid  \mathbf{Y}}\right)\) ,但从 (3.28) 来看很难直接抽样。为此,使用标准化的重要抽样法。为了解除 \(\left( {K,\eta }\right)\) 的取值限制,作变换 \(\alpha  = \log K,\beta  = \log \frac{\eta }{1 - \eta }\) ,则 \(\alpha ,\beta  \in  \left( {-\infty ,\infty }\right)\) ,而(3.28)对应的 \(\left( {\alpha ,\beta }\right)\) 的后验密度为:

\[
p\left( {\alpha ,\beta  \mid  \mathbf{Y}}\right)  \propto  \frac{{e}^{\alpha }}{{\left( 1 + {e}^{\alpha }\right) }^{2}}\mathop{\prod }\limits_{{j = 1}}^{n}\frac{B\left( {\frac{{e}^{\alpha }}{1 + {e}^{-\beta }} + {y}_{j},\frac{{e}^{\alpha }}{1 + {e}^{\beta }} + {n}_{j} - {y}_{j}}\right) }{B\left( {\frac{{e}^{\alpha }}{1 + {e}^{-\beta }},\frac{{e}^{\alpha }}{1 + {e}^{\beta }}}\right) }. \tag{3.29}
\]

取值无限制的随机变量试抽样密度经常使用自由度较小的 \(\mathrm{t}\) 分布,比如 \(\mathrm{t}\left( 4\right)\) 分布,设 \(\mathrm{t}\left( 4\right)\) 分布密度函数为 \(g\left( \cdot \right)\) ,用独立的 \(\mathrm{t}\left( 4\right)\) 分布生成 \(\left( {\alpha ,\beta }\right)\) 的试抽样样本 \(\left( {{\alpha }_{i},{\beta }_{i}}\right) ,i = 1,2,\ldots ,N\) ,可以估计 \(E\left( {\log K \mid  \mathbf{Y}}\right)\) 为

\[
\widehat{\alpha } = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{\alpha }_{i}\frac{p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right) }{g\left( {\alpha }_{i}\right) g\left( {\beta }_{i}\right) }}{\mathop{\sum }\limits_{{i = 1}}^{N}\frac{p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right) }{g\left( {\alpha }_{i}\right) g\left( {\beta }_{i}\right) }}. \tag{3.30}
\]

其中的 \(p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right)\) 只要用(3.29)的右侧计算,因为分子和分母的归一化常数可以消掉。

